library(fixest) # to demean variables
library(multiwayvcov) # to cluster SEs
library(sandwich) # to get HC2 SEs
library(stargazer) # to generate tables
library(mvtnorm) # to simulate betas and vcov matrix
library(betareg) # to run beta regression

rm(list = ls())

# load data
load(file = 'chData.RData')

# Put each dataset in a different object
ch_reform <- ch_list[[1]]
ch_refTerm <- ch_list[[2]]
ch_pre_post <- ch_list[[3]]
ch_refTermMChanges <- ch_list[[4]]
ch_refTermParty <- ch_list[[5]]
ch_refTermMChangesParty <- ch_list[[6]]

# Exclude "Gonzalo Fuenzalida" and 'Javier Macaya' because their district changed after the reform
ch_reform <- subset(ch_reform, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTerm <- subset(ch_refTerm, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_pre_post <- subset(ch_pre_post, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTermMChanges <- subset(ch_refTermMChanges, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTermParty <- subset(ch_refTermParty, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTermMChangesParty <- subset(ch_refTermMChangesParty, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))

# Generate DVs
ch_refTerm$women_pct_beta <- ch_refTerm$women_pct/100
ch_refTerm$women_pct_beta <- (ch_refTerm$women_pct_beta*(nrow(ch_refTerm) - 1) + 0.5)/nrow(ch_refTerm)

ch_refTermMChanges$women_pct_beta <- ch_refTermMChanges$women_pct/100
ch_refTermMChanges$women_pct_beta <- (ch_refTermMChanges$women_pct_beta*(nrow(ch_refTermMChanges) - 1) + 0.5)/nrow(ch_refTermMChanges)

ch_reform$women_pct_beta <- ch_reform$women_pct/100
ch_reform$women_pct_beta <- (ch_reform$women_pct_beta*(nrow(ch_reform) - 1) + 0.5)/nrow(ch_reform)

ch_pre_post$women_pct_beta <- ch_pre_post$women_pct/100
ch_pre_post$women_pct_beta <- (ch_pre_post$women_pct_beta*(nrow(ch_pre_post) - 1) + 0.5)/nrow(ch_pre_post)

# M1: G*M C1,C2,C3 (all data)
m1 <- betareg(women_pct_beta ~ magAnalysis + female + magAnalysisWoman + session, data = ch_refTerm)
vcov1 <- vcovCL(m1, cluster = ch_refTerm$bill_author)
s1  <- sqrt(diag(vcov1))

# M2: G*M C1,C2,C3 (M changes at the middle of C2)
m2 <- betareg(women_pct_beta ~ magAnalysis + female + magAnalysisWoman + session + reform, data = ch_refTermMChanges)
vcov2 <- vcovCL(m2, cluster = ch_refTermMChanges$bill_author)
s2  <- sqrt(diag(vcov2))

# M3: G*M C1,C3
small <- subset(ch_refTerm, session %in% c('2010-2014', '2018-2022'))
small$women_pct_beta <- ((small$women_pct/100)*(nrow(small) - 1) + 0.5)/nrow(small)
m3 <- betareg(women_pct_beta ~ magAnalysis + female + magAnalysisWoman + session, data = small)
vcov3 <- vcovCL(m3, cluster = small$bill_author)
s3  <- sqrt(diag(vcov3))

# M4: C2
inSample <- names(table(ch_pre_post$bill_author)[table(ch_pre_post$bill_author) == 2])
ch_pre_post <- subset(ch_pre_post, bill_author %in% inSample)
m4 <- betareg(women_pct_beta ~ magAnalysis + magAnalysisWoman + bill_author, data = ch_pre_post)
vcov4 <- vcovCL(m4, cluster = ch_pre_post$bill_author)
s4  <- sqrt(diag(vcov4))

# M5: G*M C1,C2,C3 (only repeated legislators)
small <- subset(ch_refTerm, allSessions == 1)
small$women_pct_beta <- ((small$women_pct/100)*(nrow(small) - 1) + 0.5)/nrow(small)
m5 <- betareg(women_pct_beta ~ magAnalysis + magAnalysisWoman + bill_author + session, data = small)
vcov5 <- vcovCL(m5, cluster = small$bill_author)
s5  <- sqrt(diag(vcov5))

# M6: G*M C1, C2, C3 (M changes at the middle of C2) (only repeated legislators)
small <- subset(ch_refTermMChanges, allSessions == 1)
small$women_pct_beta <- ((small$women_pct/100)*(nrow(small) - 1) + 0.5)/nrow(small)
m6 <- betareg(women_pct_beta ~ magAnalysis + magAnalysisWoman + session + reform + bill_author, data = small)
vcov6 <- vcovCL(m6, cluster = small$bill_author)
s6  <- sqrt(diag(vcov6))

# M7: G*M C1,C3 (only repeated legislators)
small <- subset(ch_refTerm, session %in% c('2010-2014', '2018-2022') & l10_l18 == 1)
small$women_pct_beta <- ((small$women_pct/100)*(nrow(small) - 1) + 0.5)/nrow(small)
m7 <- betareg(women_pct_beta ~ magAnalysis + magAnalysisWoman + bill_author + session, data = small)
vcov7 <- vcovCL(m7, cluster = small$bill_author)
s7  <- sqrt(diag(vcov7))

# M8: C3
small <- subset(ch_refTerm, session == '2018-2022')
small$women_pct_beta <- ((small$women_pct/100)*(nrow(small) - 1) + 0.5)/nrow(small)
m8 <- betareg(women_pct_beta ~ magAnalysis + female + magAnalysisWoman, data = small)
vcov8 <- vcov(m8)
s8  <- sqrt(diag(vcov8))

# Table F.27
stargazer(m1, m2, m3, m4, m5, m6, m7, m8, 
	se = list(s1, s2, s3, s4, s5, s6, s7, s8), 
	covariate.labels = c('M', 'Woman', 'M x Woman'), 
	add.lines = list(c('FE by Legislative Term', 'Yes', 'Yes', 'Yes', 'No', 'Yes', 'Yes', 'Yes', 'No'),
		c('FE by Reform', 'No', 'Yes', 'No', 'No', 'No', 'Yes', 'No', 'No'),
		c('FE by Legislator', 'No', 'No', 'No', 'Yes', 'Yes', 'Yes', 'Yes', 'No')),
	dep.var.labels = "Bills' Portfolio on Women's Issues (%)",
	omit.stat = c("ser",'f'),
	omit = c('reform', 'bill_author', 'session'),
	star.cutoffs = c(0.10, 0.05, 0.01),
	no.space = TRUE,
	single.row = FALSE,	
	notes.append = FALSE,
	title = "Association Between Legislative Portfolio, Gender, and District Magnitude--2014-2021  Cámara de Diputadas y Diputados de Chile",
	label = 't:mag'
)

#################### Substantive Effect ########################
## Scenarios: Women, M = 2
W2 <- data.frame(1, mag = 2, female = 1, magAnalysisWoman = 2 * 1)
## Scenarios: Women, M = 3 (post-reform)
W3 <- data.frame(1, mag = 3, female = 1, magAnalysisWoman = 3 * 1)
## Scenarios: Men, M = 2
M2 <- data.frame(1, mag = 2, female = 0, magAnalysisWoman = 2 * 0)
## Scenarios: Men, M = 3 (post-reform)
M3 <- data.frame(1, mag = 3, female = 0, magAnalysisWoman = 3 * 0)
# Women and Men Data
Wdata <- rbind(W2, W3)
Mdata <- rbind(M2, M3)

## Sims
set.seed(123)
sims1 <- rmvnorm(n = 1000, mean = coef(m1), sigma = vcov1)[,c('(Intercept)', 'magAnalysis', 'female', 'magAnalysisWoman')]
sims2 <- rmvnorm(n = 1000, mean = coef(m2), sigma = vcov2)[,c('(Intercept)', 'magAnalysis', 'female', 'magAnalysisWoman')]
sims3 <- rmvnorm(n = 1000, mean = coef(m3), sigma = vcov3)[,c('(Intercept)', 'magAnalysis', 'female', 'magAnalysisWoman')]
sims4 <- rmvnorm(n = 1000, mean = coef(m4), sigma = vcov4)[,c('(Intercept)', 'magAnalysis', 'magAnalysisWoman')]
sims5 <- rmvnorm(n = 1000, mean = coef(m5), sigma = vcov5)[,c('(Intercept)', 'magAnalysis', 'magAnalysisWoman')]
sims6 <- rmvnorm(n = 1000, mean = coef(m6), sigma = vcov6)[,c('(Intercept)', 'magAnalysis', 'magAnalysisWoman')]
sims7 <- rmvnorm(n = 1000, mean = coef(m7), sigma = vcov7)[,c('(Intercept)', 'magAnalysis', 'magAnalysisWoman')]
sims8 <- rmvnorm(n = 1000, mean = coef(m8), sigma = vcov8)[,c('(Intercept)', 'magAnalysis', 'female', 'magAnalysisWoman')]


## Calculate Predicted Values
# Women
pred1_W <- apply(plogis(sims1 %*% t(Wdata[2,])) - plogis(sims1 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025)) * 100
pred2_W <- apply(plogis(sims2 %*% t(Wdata[2,])) - plogis(sims2 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025)) * 100
pred3_W <- apply(plogis(sims3 %*% t(Wdata[2,])) - plogis(sims3 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025)) * 100
pred4_W <- apply(plogis(sims4 %*% t(Wdata[2, -3])) - plogis(sims4 %*% t(Wdata[1, -3])), 2, quantile, c(0.975, 0.5, 0.025)) * 100
pred5_W <- apply(plogis(sims5 %*% t(Wdata[2, -3])) - plogis(sims5 %*% t(Wdata[1, -3])), 2, quantile, c(0.975, 0.5, 0.025)) * 100
pred6_W <- apply(plogis(sims6 %*% t(Wdata[2, -3])) - plogis(sims6 %*% t(Wdata[1, -3])), 2, quantile, c(0.975, 0.5, 0.025)) * 100
pred7_W <- apply(plogis(sims7 %*% t(Wdata[2, -3])) - plogis(sims7 %*% t(Wdata[1, -3])), 2, quantile, c(0.975, 0.5, 0.025)) * 100
pred8_W <- apply(plogis(sims8 %*% t(Wdata[2,])) - plogis(sims8 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025)) * 100

# Men
pred1_M <- apply(plogis(sims1 %*% t(Mdata[2,])) - plogis(sims1 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025)) * 100
pred2_M <- apply(plogis(sims2 %*% t(Mdata[2,])) - plogis(sims2 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025)) * 100
pred3_M <- apply(plogis(sims3 %*% t(Mdata[2,])) - plogis(sims3 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025)) * 100
pred4_M <- apply(plogis(sims4 %*% t(Mdata[2, -3])) - plogis(sims4 %*% t(Mdata[1, -3])), 2, quantile, c(0.975, 0.5, 0.025)) * 100
pred5_M <- apply(plogis(sims5 %*% t(Mdata[2, -3])) - plogis(sims5 %*% t(Mdata[1, -3])), 2, quantile, c(0.975, 0.5, 0.025)) * 100
pred6_M <- apply(plogis(sims6 %*% t(Mdata[2, -3])) - plogis(sims6 %*% t(Mdata[1, -3])), 2, quantile, c(0.975, 0.5, 0.025)) * 100
pred7_M <- apply(plogis(sims7 %*% t(Mdata[2, -3])) - plogis(sims7 %*% t(Mdata[1, -3])), 2, quantile, c(0.975, 0.5, 0.025)) * 100
pred8_M <- apply(plogis(sims8 %*% t(Mdata[2,])) - plogis(sims8 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))  * 100


# Building blocks of table 2
bill_hat_women <- c(round(pred1_W[2,], 2), round(pred2_W[2,], 2), round(pred3_W[2,], 2),
	round(pred4_W[2,], 2), round(pred5_W[2,], 2), round(pred6_W[2,], 2),
	round(pred7_W[2,], 2), round(pred8_W[2,], 2))
ci_women <- c(paste0('(', round(pred1_W[3,], 2), ', ', round(pred1_W[1,], 2), ')'), 
	paste0('(', round(pred2_W[3,], 2), ', ', round(pred2_W[1,], 2), ')'), 
	paste0('(', round(pred3_W[3,], 2), ', ', round(pred3_W[1,], 2), ')'), 
	paste0('(', round(pred4_W[3,], 2), ', ', round(pred4_W[1,], 2), ')'), 
	paste0('(', round(pred5_W[3,], 2), ', ', round(pred5_W[1,], 2), ')'), 
	paste0('(', round(pred6_W[3,], 2), ', ', round(pred6_W[1,], 2), ')'), 
	paste0('(', round(pred7_W[3,], 2), ', ', round(pred7_W[1,], 2), ')'), 
	paste0('(', round(pred8_W[3,], 2), ', ', round(pred8_W[1,], 2), ')'))


bill_hat_men <- c(round(pred1_M[2,], 2), round(pred2_M[2,], 2), round(pred3_M[2,], 2),
	round(pred4_M[2,], 2), round(pred5_M[2,], 2), round(pred6_M[2,], 2),
	round(pred7_M[2,], 2), round(pred8_M[2,], 2))
ci_men <- c(paste0('(', round(pred1_M[3,], 2), ', ', round(pred1_M[1,], 2), ')'), 
	paste0('(', round(pred2_M[3,], 2), ', ', round(pred2_M[1,], 2), ')'), 
	paste0('(', round(pred3_M[3,], 2), ', ', round(pred3_M[1,], 2), ')'), 
	paste0('(', round(pred4_M[3,], 2), ', ', round(pred4_M[1,], 2), ')'), 
	paste0('(', round(pred5_M[3,], 2), ', ', round(pred5_M[1,], 2), ')'), 
	paste0('(', round(pred6_M[3,], 2), ', ', round(pred6_M[1,], 2), ')'), 
	paste0('(', round(pred7_M[3,], 2), ', ', round(pred7_M[1,], 2), ')'), 
	paste0('(', round(pred8_M[3,], 2), ', ', round(pred8_M[1,], 2), ')'))

tableF28 <- rbind("Women "= bill_hat_women, " " = ci_women, 
	"Men "= bill_hat_men, " " = ci_men)
colnames(tableF28) <- paste0('Model ', 1:8)
stargazer(tableF28,
	title = "Predicted percentage point change in the cosponsorship portfolio on women’s issues when M increases by one")

